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1.  Introduction 


Dielectric  materials  with  charged  defects  exhibit  a  variety  of  physical  phenomena  whose  origins 
are  not  fully  understood.  Early  continuum  theories  (Devonshire,  1954;  Toupin,  1956)  of  the 
electromechanical  behavior  of  dielectric  media  have  been  set  forth,  though  these  do  not  consider 
defects  explicitly.  However,  defects  such  as  vacancies  have  been  the  focus  of  considerable  study 
(Lifshitz,  1963),  particularly  with  regards  to  crystalline  ceramics  of  current  interest  to  the  U.S. 
Army  Research  Laboratory. 

The  unique  aspect  of  the  present  study  is  consideration  of  electrically  charged,  as  opposed  to 
neutral,  vacancies  in  dielectric  solids.  Defect  concentrations  and  excess  charges  in  ferroic 
ceramics  may  be  adjusted  during  processing  via  heat  treatments  and/or  addition  of  doping 
chemicals  (Cole  et  ah,  2003).  Charged  point  defects  have  been  identified  as  a  major  factor 
affecting  the  reliability  of  ferroelectric  devices  (Damjanovic,  1998),  including  gate  dielectric 
semiconductors,  particularly  those  of  thin  film  geometry  (Buchanan,  1999). 

Here,  a  general  modeling  framework  is  constructed  for  elastic  dielectric  semiconductors  with 
mobile  charged  point  vacancies.  This  framework  combines  the  physics  of  continuum  elasticity, 
electrostatics,  mass  diffusion,  and  charged  defect  kinetics.  Changes  in  surface  morphology  due 
to  the  boundary  flux  of  charged  vacancies  are  captured,  extending  a  previous  theory  of  one  of  the 
co-authors  on  neutral  vacancy  kinetics  (Grinfeld  and  Hazzledine,  1997). 

The  theory  is  implemented  numerically  in  a  finite  difference  code  (Hoffman,  1992)  enabling 
simultaneous  solution  of  the  elliptic  equations  of  electrostatics  of  dielectrics  and  the  transient 
parabolic  equations  of  charged  diffusion.  The  analysis  is  limited  to  a  single  spatial  dimension. 
The  time  duration  of  the  problem  is  decomposed  into  a  sequence  of  steps.  A  second-order 
accurate  fully  implicit  scheme  is  invoked  to  solve  Maxwell's  equations  in  each  step,  while  a  fully 
explicit  scheme  is  used  to  integrate  the  transient  vacancy  concentration.  The  spatial  domain  and 
grid  spacing  are  updated  when  the  surface  flux  of  concentration  is  nonzero,  as  vacancies  exiting 
the  domain  influence  its  instantaneous  dimensions. 

Documentation  is  presented  for  the  computer  implementation.  Included  here  are  descriptions  of 
the  source  code  structure,  user  instructions,  and  representative  input  fdes  for  the  software,  the 
latter  specifically  for  analysis  of  barium  strontium  titanate  (Bai.xSrxTiOs)  (BST)  thin  films 
containing  charged  oxygen  vacancies.  The  source  code  is  given  in  the  appendix. 

In  the  notation  that  follows,  the  Einstein  summation  convention  is  used  on  repeated  lower-case 
indices,  unless  indicated  otherwise.  Cartesian  spatial  coordinate  indices  span  three  dimensions 
and  are  written  in  Roman  font,  while  curvilinear  surface  coordinate  indices  span  two  dimensions 
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and  are  written  in  Greek  font.  Subseripted  eommas  denote  eovariant  differentiation.  Capitalized 
subseripts  are  often  used  in  physical  constants  and  are  not  summed.  In  the  description  of 
numerical  methods,  subscripts  are  often  used  for  node  numbers,  and  superscripts  for  time 
increments. 


2.  Theory 


The  governing  equations,  thermodynamic  framework,  and  constitutive  relations  are  first 
presented,  with  limited  derivations,  in  three-dimensional  (3-D)  form  in  section  2.1.  The  one¬ 
dimensional  (1-D)  equations  that  are  solved  numerically  follow  in  section  2.2. 

2,1  Model  Framework 

The  local  electrostatic  behavior  of  dielectric  continua  is  dictated  by  Maxwell's  equations 
(Stratton,  1941): 


D)=P, 

(1) 

1 

II 

(2) 

(3) 

where  D  is  the  electric  displacement,  p  is  the  charge  density,  E  is  the  electric  field,  (j)  is  the 
electrostatic  potential,  and  is  the  permittivity  of  free  space.  The  polarization  vector  P  is 

defined  only  within  the  material  and  vanishes  in  free  space.  Local  mechanical  equilibrium  and 
mass  conservation  are  ensured  by 

r75=0  (4) 

and 

p  +  pw)  =0,  (5) 

where  <j  is  the  symmetric  Cauchy  stress,  p  is  the  mass  density  of  the  deformed  solid  with 
vacancies,  and  u  is  the  displacement,  which  itself  encompasses  expansion  or  contraction  due  to 
vacancies  within  the  bulk  material.  Small  displacements  are  assumed  henceforth  in  the  present 
theory.  The  balance  of  energy  and  the  dissipation  inequality  may  be  written  in  global  rate  form 
as 

0  =  w-n>o,  (6) 

with  0  the  dissipation,  W  the  rate  of  external  work,  and  E2  the  system  energy.  For  a  dielectric 
solid  containing  mobile  charged  vacancies,  the  external  work  and  energy  of  the  system  are 
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(V) 


and 


^  =  —j  ^<i)ds  +  J t'ii-ds  - 1  juQ'n-ds 


dt 


i?  =  J  p{u  -6i)^dv 


(8) 


where  s  and  v  are  the  surfaee  and  volume  of  the  system  with  unit  normal  n,  a  is  the  surfaee 
eharge  density,  i  is  the  traetion  veetor,  ju  is  the  ehemieal  potential  for  vaeaney  diffusion,  Q  is 
the  vaeaney  flux,  U  is  the  internal  energy  per  unit  mass  of  the  substanee,  and  7  is  the  surfaee 
energy  density.  Heat  eonduetion  is  not  eonsidered  explieitly  here,  as  is  evident  from  equations 
6-8.  The  Helmholtz  free  energy  y/ ,  speeifie  entropy  77 ,  and  temperature  6  are  related  by  the 
usual  thermo  dynamie  relationship 

xj/  =  U  -rjO .  (9) 


The  following  eonstitutive  assumptions  are  made  regarding  the  free  energy  and  eharge  density: 


y/  =  \l/(u^  .,P' ,6,^'], 


(10) 


and 


p  =  ez^ , 


(11) 


where  ^  is  the  number  of  vaeaneies  per  unit  volume.  Assumption  (equation  10)  suggests  a  free 
energy  dependenee  on  meohanieal  strain,  polarization,  temperature,  and  vaeaney  eoneentration. 
Equation  1 1  denotes  that  the  eharge  density  is  proportionate  to  the  vaeaney  concentration,  with  e 
and  z  the  charge  of  an  electron  and  the  valence  contribution  of  each  defect,  respectively. 
Substituting  equations  9-1 1  into  equation  6  and  making  use  of  the  balance  laws  in  equations  1-5 
and  the  divergence  theorem,  and  restricting  ^  =  0  on  5,  the  following  thermodynamically 

admissible  bulk  constitutive  relations  may  be  derived: 


E.  =  pdy/ jdP' , 

(12) 

cr"  ^  p By/ jdu^^  , 

(13) 

rj  =  -dy/jdO , 

(14) 

(15) 

Q' 

(16) 

p  =  pdyr  jd^  ez(j) , 

(17) 
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where  the  diffusivity  tensor  d  is  symmetrie  and  positive  definite.  Let  j  denote  the  volume 
fraction  of  vacancies,  with  «  is  a  scalar  conversion  factor.  Mass  conservation  requires  that 

=  (18) 

where  d  is  the  surface  velocity  of  the  substance  with  unit  normal  n  and  q  is  the  flux  tangential  to 
the  surface  of  charged  vacancies.  Conservation  of  mechanical  work  leads  to 

i' =  pdy/ |^u^  jHj .  (19) 

By  requiring  that  the  volumetric  energy  of  the  system  O  remain  constant  or  decrease  with  time 
in  the  absence  of  mechanical  working  due  to  external  stresses,  the  following  relations  for  the 
in-plane  and  normal  surface  fluxes  are  formulated: 

q-=A-^({p^  +  p^x)l(l-x))j,.  (20) 

and 

[\- x)Q'n.=  pi^p\l/  + p(l)  +  [\- x)[pdy/ldx  +  a'ez(l)fj,  (21) 

where  A  is  the  symmetric  positive  definite  surface  diffusivity  tensor  and  is  a  positive  scalar 
characterizing  the  resistance  of  the  surface  to  penetration  by  vacancies.  Implicit  in  equations  4, 
13,  and  19  is  the  vanishing  of  Maxwell's  stress  tensor  (Toupin,  1956),  meaning  that  contributions 
from  terms  of  second  order  in  the  electric  field  and  polarization  are  neglected  in  the  mechanical 
equilibrium  equations. 

The  continuum  theory  is  applied  here  to  describe  linear  elastic  dielectric  solids.  A  specific  form 
of  the  free  energy  density  is  thus  postulated  as 

=(1/2)C»%  +(II2)X,FP‘  +<p(i,0) ,  (22) 

where  C  and  X  are  linear-elastic  moduli  and  inverse  dielectric  susceptibility,  respectively,  and 
symmetrized  indices  are  in  parentheses.  Presently,  we  address  only  the  response  of  the  material 
in  its  paraelectric  state,  at  temperatures  above  the  Curie  point.  Thus,  phase  transformations, 
piezoelectricity,  pyroelectricity,  and  spontaneous  polarization  are  not  considered.  From  equation 
22,  bulk  thermodynamic  relations  in  equations  12,  13,  and  3  reduce  to 


(23) 

(24) 

II 

(25) 
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with  £■!  =  d‘^  +£’o'/l the  relative  permittivity  (i.e.,  real  dielectrie  constant).  The  vacancy  and 
temperature-dependent  contribution  (p  is  assumed  here  to  follow  the  universal  relation  for  the 
chemical  potential  of  an  ideal  mixture  (Fried  et  ah,  1977),  most  relevant  for  noninteracting 
species  and  small  vacancy  concentrations: 

(pN A  (0)  +  NAk^Ov  In  v) .  (26) 

In  equation  26,  Na  is  Avagadro’s  number,  kg  is  Boltzmann’s  constant,  Nj.=\/a  is  the  number 
of  atomic  sites  per  unit  volume,  v  is  the  mole  fraction  of  vacancies,  and  Gq  is  the  bulk  Gibbs  free 
energy  of  unstressed,  defect-free  substance.  As  implemented  here,  we  assume  simply  that 

G,  =  e[c{e)-ri[e)\,  (27) 

with  specific  heat  c  and  specific  entropy  fj .  Upon  application  of  equation  26,  the  kinetic 
equation  for  bulk  diffusion  then  follows  directly  from  equation  16  as 

a=-d‘^[kgexJx-ezE^).  (28) 


2,2  One-Dimensional  Model 

The  physical  system  analyzed  is  a  substance  of  length  (i.e.,  thickness)  T  with  applied  bias 
voltage  V,  as  shown  in  figure  1 .  Alternatively  to  voltage  boundary  conditions,  electric  fields  may 
be  applied  at  the  boundaries  (not  shown  in  figure  1).  No  mechanical  tractions  are  applied,  and 
isotropic  material  properties  are  assumed  {s'g=  Sgd'^ , =  dS'^ ),  thus  rendering  the  analysis  one 
dimensional  in  thickness  direction  x.  Correspondingly  in  this  analysis,  surface  flux  q“  =0  and 
surface  tension  r  =  0 .  A  constant  temperature  0  is  assumed. 


Under  the  preceding  conditions,  the  governing  electrostatic  equations  1-3  reduce  to 

dDjdx  =  Sf^SgdE  ldx=-  s^Sg  {d^^l  dx^ P  ■  (29) 

The  vacancy  evolution  equations  28  and  15  become 
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and 


Q  =  -d(  kg0 1  z{dz/dx)-ezE) 


(30) 


Z  =  -a{dQldx),  (31) 

with  the  domain  size  evolving,  from  equation  18,  as 

a  =  a{Q(T)l(\- x(T))-Q(0)l(\- x(f>)))-  (32) 

Spatial  boundary  conditions  are  applied  as  follows: 

i{0,t)  =  i{L,t)  =  0,  (33) 

^(0,t)  =  Fo  or  E{Q,t)  =  E„  (34) 

(p{T,t)  =  V^  or  E{T,t)  =  E^,  (35) 

x{^d)  =  Xo  or  0(O,t)  =  go.  (36) 

and 

x{T,t)  =  XT  or  Q{T,t)  =  QT-  (37) 

The  boundary  fluxes  gg  and  Qj.  may  be  prescribed  as  constant  values,  or  determined 
instantaneously  from  the  model  thermodynamics  (equation  21): 

=  +  +  +  «  W))-  (38) 

Initial  conditions  for  the  vacancy  concentration  and  electrostatic  charge  density  are  applied  as 

xix,0)  =  z{x)  (39) 

and 

p{x,0)  =  p{x),  (40) 

where  p  =  ezz  /  cc  should  be  imposed  for  consistency  with  equation  1 1 . 

Note  from  equations  30  and  31  that  the  electric  field  may  act  as  a  local  source  or  sink  for 
vacancies,  implying  that  ions  may  be  impinged  within  or  released  from  the  lattice  due  to 
electrostatic  forces.  In  some  cases,  it  is  convenient  and  realistic  to  impose  an  additional 
constraint  that  the  total  number  of  vacancies  in  the  system  remains  constant,  i.e., 

f^Xdx^Zj,  (41) 

where  Zo  the  initial  average  concentration. 
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3.  Numerical  Methods 


The  time  duration  of  the  problem  is  deeomposed  into  a  series  of  steps.  For  each  time  step,  the 
electrostatic  problem  is  solved,  with  the  corresponding  solution  for  the  electric  field  used  in  the 
transient  equations  (equation  30)  for  updating  the  vacancy  concentration.  In  the  following  time 
step,  this  updated  vacancy  concentration  is  used  to  compute  the  local  electric  charge  entering 
Maxwell’s  equations  (equation  29).  The  analysis  thus  marches  forward  in  time,  with  the 
numerical  solutions  of  the  equations  of  electrostatics  and  diffusion  coupled  in  this  manner. 

The  1-D  spatial  domain  is  discretized  into  n-\  increments  of  equal  length  Ax ,  with  n  the  total 
number  of  nodes,  as  shown  in  figure  2. 


3.1  Electrostatics 

A  second-order  accurate  implicit  scheme  is  used  to  solve  Maxwell's  static  equations  in  each  step. 
Within  the  domain  0  <  x  <  T  ,  equation  29  is  written,  to  accuracy  of  order  (Ax)^ ,  using  a  finite 
difference  approximation  (Hoffman,  1992)  for  the  second  spatial  derivative  of  (j) ; 

-  2^,  +  =  -Pi  (Ax)'  /  ,  (42) 

where  subscripts  denote  node  numbers.  The  solution  of  equation  42  for  all  nodes  i  is  then 
written  in  matrix  form  as 

M  =  [A]->],  (43) 

where  [^]  is  the  n-dimensional  solution  vector  of  nodal  electric  potentials,  [  A]  is  the  n  x  n 
coefficient  matrix  that  is  generally  tridiagonal  and  sparse,  and  [/?]  is  an  n-dimensional  vector 

that  corresponds  to  the  right  hand  side  of  equation  42.  In  the  solution  procedure,  the  dimensions 
of  [<!>],  [A]  ,  and  [ /?]  are  reduced  by  one  for  each  electrostatic  potential  boundary  condition 

enforced  a  priori.  Once  the  electric  potential  is  known,  the  electric  field  is  then  computed 
numerically  as 
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(44) 


The  following  second-order  finite  difference  approximations  are  used  at  the  boundaries  (/  =  1 
and  i  =  n)\ 

-llSxE^  =  -  <t>-i ,  (45) 

and 

-lExE^  =  (t>n-l  -  +  '^<t>n  '  (46) 

3.2  Transient  Diffusion 

A  fully  explicit  scheme  is  used  to  integrate  the  transient  vacancy  concentration,  with  all 
necessary  spatial  derivatives  obtained  using  a  second-order  accurate  finite  difference  approach. 
The  vacancy  flux  (equation  30)  is  written  as  follows,  where  subscripts  denote  node  numbers; 

Q^=-d[d[pd  y/ 1  dXi )!  dx-  ezE^ ) ,  (47) 

where  for  nodes  in  the  domain  0  <  x  <  T  , 

pd\j/ldx,=k^e\x^Xi  (48) 

and 

d{pdy/ldXi)l dx  =  kge{\vL  -  in  ) / 2 Ax  .  (49) 

Note  that  in  equation  47,  the  electric  field  E.  is  obtained  from  the  electrostatic  solution  given  by 

equations  44-46.  At  the  boundaries,  the  following  approximations  were  used,  analogous  to 
equations  45  and  46: 

d{pdy//dXi)/dx  =  kg6{-{nx^+^^r).X2  -31n  jj)/2Ax,  (50) 

and 

flf(/75^^/aj„)/ix  =  ^^^(lnj„_2-41nj„_i+31njJ/2Ax.  (51) 

The  rate  equation  for  concentration  is  approximated  as  follows  for  0<x<T  : 

=  -« ( Qm  -  Qi  i )  /  2Ax ,  (52) 

and  at  the  boundaries  as 

Xi=-a{-Q^+4Q2-3Qi)/2Ax  (53) 

and 

=-«(a-2-40„_i+30„)/2Ax.  (54) 

Application  of  concentration  boundary  conditions  simply  entails  =  0  ,  while  application  of 
constant  flux  conditions  at  a  boundary  eliminates  the  need  to  solve  (equation  47)  at  that 
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boundary.  When  transient  flux  eonditions  (equation  38)  are  applied,  the  following 
approximation  suffiees: 

(1  -  Z )  Qi  =  -P [P¥i  +  pA  +  (1  -  Z  )  ksO In  ) ,  (55) 

where  the  free  energy  per  unit  volume  is  found  from 

P¥i  =  \Ei I V 2  +  (Go  +  N^kge{Xi  > ) In (z  > ))I^a-  (56) 

Finally,  the  eoneentration  is  updated  explieitly  as 

zr  =Z;+ZAt,  (57) 


where  At  is  the  time  inerement. 


The  logarithmie  form  of  the  chemieal  potential,  and  flux  equations  48-51,  require  that  >  0  . 
This  is  achieved  in  practice  by  enforcing 

(58) 


where  is  projected  updated  concentration  from  equation  57  and  Xmm  -  iO  is  n  near- 
negligible,  default  minimum  concentration.  Optional  global  conservation  condition  (equation 
41)  is  imposed  by  rescaling  equation  57; 


Z, 


n-\ 


=  XqT [x‘i  +  z  At)  - +  Z  At) 


(59) 


The  domain  size  is  updated  from  equation  32  as 

(60) 


and 


Ax‘^^‘  ={T  +  dAt)/{n-l). 


(61) 


For  the  solution  of  equation  42  in  the  next  time  step,  the  updated  charge  concentration  is  found 
from  equation  1 1 : 

Pi  =ezXi  la.  (62) 

A  fixed  time  increment  is  used  throughout  the  analysis,  i.e., 

=  (63) 

where  and  m  are  the  total  simulation  time  and  number  of  time  increments,  respectively. 
Convergence  and  stability  of  the  numerical  solution  dictate  the  practical  choice  of  At 
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(Chapra  and  Canale,  1998).  Let  D  =  dkj^O  I  ,  where  is  a  referenee  eoneentration  per  unit 
volume,  on  the  order  of  the  minimum  loeal  eoneentration  in  the  domain.  Then  one  may  write 

M<A{^xflD,  (64) 

where  A<\I2  guarantees  convergenee  and  stability.  Setting  A<\l  A  ensures  that  the  solution 
will  not  oseillate,  while  the  ehoiee  A  =  \l  6  has  been  shown,  for  the  pure  diffusion  problem,  to 
minimize  truneation  error  (Carnahan  et  ah,  1969).  A  deseription  of  what  effeet  the  unique 
eleetric  field  term  ezE.  in  equation  47  may  have  on  stability  and  eonvergenee  is  not  presently 

available. 


4.  Software  Manual 


In  what  follows,  the  souree  eode  strueture,  input  fdes,  and  output  fdes  are  doeumented.  Note 
that  for  referenee,  the  complete  code  is  contained  in  the  appendix. 

4,1  Code  Structure 

The  source  code  is  written  in  the  standard  FORTRAN90  language.  It  has  been  compiled  with 
Portland  Group's  FORTRAN  compiler  (i.e.,  pgf90)  (http://www.pgroup.com/)  and  executed  on  a 
64-bit  Linux  workstation.  Floating  point  values  are  represented  in  scientific  double  precision. 
This  is  a  serial,  as  opposed  to  parallel,  code.  The  program  permits  solution  of  the  elliptic 
differential  equations  of  electrostatics  and/or  the  parabolic  differential  equations  for  transient 
diffusion.  Standard  SI  units  are  used  throughout.  The  software  enables  all  of  the  features 
described  in  section  3,  as  well  as  a  few  additional  options  offering  additional  user  flexibility. 

The  code  consists  of  the  following  routines: 

•  the  main  routine,  'semiconductor  ! D',  which  controls  the  grid  spacing,  memory  allocation, 
and  time  incrementation. 

•  subroutine  'elliptic',  which  solves  Maxwell's  equations. 

•  subroutine  'parabolic',  which  solves  the  transient  diffusion  equations. 

•  subroutine  'chemjotential',  which  computes  the  local  free  energy  and  its  gradient. 

•  subroutine  'LU_decomp',  which  decomposes  the  [A]  matrix  in  (43)  to  LU-form. 

•  subroutine  'LU  backsub',  which  solves  equation  43  using  a  back  substitution  technique. 
Figure  3  is  a  flowchart  demonstrating  execution  of  the  code. 
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semiconduclor_1D 
-read  'parameters.inp’ 

-loop  for  t<tm 

-if  electrostatic  problem,  call  ‘elliptic’ 
-if  diffusion  problem,  call  ‘parabolic’ 
-write  ‘carray.bd’  and  ‘earray.bd’ 

-end  loop 


elliptic 

parabolic 

-read  ‘electrostatic.inp’ 

-update  X 

-update  charge 

-read  ‘diffusion.inp’ 

-compute  [A],  [p] 

-call  ‘chem_potentiar 

-call  ‘LU_decomp’ 

-compute  0 

-call  ‘LU-backsub’ 

-compute  X,  a 

-obtain  M,  £ 

-write  ‘outdiff.bct’ 

-write  ‘outelec.brt’ 

Figure  3.  Flowchart  for  code  execution. 


4,2  Input 

Up  to  three  input  files  are  required.  The  file  read  by  the  main  routine,  'parameters.inp',  is  always 
required.  Its  format  is  as  follows: 

'parameters.inp' 

Line  1 :  problem  type.  Specify  'maxwell'  to  solve  only  the  electrostatic  equations,  'fick'  to  solve 
only  the  diffusion  equations,  or  'mixed'  to  solve  the  fully  coupled  electrostatic-diffusion 
problem 

Line  2:  number  of  nodes  n 
Line  3:  solution  end  time  [s] 

Line  4:  fixed  time  increment  At  [s] 

Line  5:  number  of  time  steps  in  between  each  write  to  the  output  files, 

Line  6:  domain  thickness  T  [m] 

Line  7:  temperature  6  [K] 

Line  8:  relative  dielectric  permittivity  ^ 

Line  9:  factor  k  enabling  permittivity  to  depend  on  local  electric  field  (Johnson,  1962): 

~  (1  +  /^^  j 

Line  10:  diffusivity  c/ 

Line  1 1 :  diffusivity  d  for  electrostatic  contribution  (set  d  =  d  ^or  consistency  with  equation 
47): 

g.  =-d{di^pdy/  ldXi)l  dx^  +  dezE. 

—3  —6  —3 

Line  12:  factor  for  conversion  of  concentration  from  [ppm]  to  [m  ],  equal  to  10  /a  [m  ] 

Line  13:  factor  for  conversion  of  concentration  from  [ppm]  to  volume  fraction,  nominally  10  ® 
Line  14:  valence  contribution  per  vacancy,  z 
Line  15:  specific  entropy  fj  [J/mol  K] 

Line  16:  specific  heat  capacity  c  [J/mol  K] 
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The  file  'electrostatie.inp'  is  read  only  if  the  problem  type  is  'maxwell'  or  'mixed'.  Its  format  is 
specified  as  follows: 

'electrostatie.inp' 

Line  1:  boundary  condition  at  x  =  0  .  Use  'potential'  to  supply  a  voltage  (j) ,  or  'flux'  to  supply 
an  electric  field  E 

Line  2:  value  of  (l>[x  =  O)  [V]  for  a  'potential'  condition,  or  value  of  £'(x  =  O)  [V/m]  for  a 
'flux'  condition 

Line  3:  boundary  condition  at  x  =  T .  Use  'potential'  to  supply  a  voltage  (j) ,  or  'flux'  to  supply 
an  electric  field  E 

Line  4:  value  of  (/)[x  =  T)  [V]  for  'potential'  condition,  or  value  of  £'(x  =  f)  [V/m]  for  'flux' 
condition 

Line  5:  text  heading  (not  used) 

Line  6:  initial  charge  density  at  node  1,  [C/m^] 

Line  7:  initial  charge  density  at  node  2,  [C/m^] 


Line  n  +  5\  initial  charge  density  at  node  n,  p\^^  [C/m^] 


The  file  'diffusion. inp'  is  read  only  if  the  problem  type  is  'fick'  or  'mixed'.  Its  format  is  given  as 
follows: 


'diffusion,  inp' 

Line  1:  boundary  condition  at  x  =  0  .  Use  'potential'  to  supply  a  constant  vacancy 
concentration  c ,  'flux'  to  supply  a  constant  vacancy  flux  Q,  or  'special'  to  apply  equation  55. 
Line  2:  value  of  c(x  =  O)  [ppm]  for  a  'potential'  condition,  value  of  Q{x  =  O)  [ppm  m/s]  for  a 
'flux'  condition,  or  value  of  [5  in  (55)  for  'special'  condition 
Line  3:  boundary  condition  at  x  =  T .  Use  'potential'  to  supply  a  constant  vacancy 
concentration  c ,  'flux'  to  supply  a  constant  vacancy  flux  Q,  or  'special'  to  apply  equation  55. 
Line  4:  value  of  c[x  =  T)  [ppm]  for  a  'potential'  condition,  value  of  Q[x  =  T)  [ppm  m/s]  for  a 
'flux'  condition,  or  value  of  [5  in  (55)  for  'special'  condition 
Line  5:  flag  to  apply  constraint  (59):  'yes'  or  'no' 

Line  6:  initial  vacancy  concentration  at  node  1,  c\Z^  [ppm] 

Line  7:  initial  vacancy  concentration  at  node  2,  c\zl  [ppm] 


Line  n  +  5\  initial  vacancy  concentration  at  node  n,  [ppm] 
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4,3  Output 


Up  to  four  output  files  are  generated,  eaeh  eontaining  eomputational  results.  The  file  'outelee.txt' 
is  written  for  problems  of  type  'maxwell'  or  'mixed',  and  provides  the  following: 


'outelee.txt' 

At  output  request  timest  =  Q,t  =  t^^J n„t  =  I n,,t  =  : 

Column  1 :  node  i 
Column  2:  time  t  [s] 

Column  3:  position  x  [m] 

Column  4:  eleetrie  potential  [V] 

Column  5:  eleetrie  field  [V/m] 

Column  6:  eleetrie  eharge  p.  [C/m^] 


The  file  'outdiff.txf  is  written  for  problems  of  type  'fiek'  or  'mixed',  and  provides  the  following: 


'outdiff.txf 

At  output  request  times  t  =  0,t  =  t^^l  n„t  =  I  =  {n,-\)t^^/  n„t  = 

Column  1 :  node  i 
Column  2:  time  t  [s] 

Column  3:  position  x  [m] 

Column  4:  vaeaney  eoneentration  c  [ppm] 

Column  5:  eoneentration  rate  c  [ppm/s] 

Column  6:  vaeaney  flux  Q  [ppm  m/s] 


The  file  'earray.txf  eontains  the  eleetrie  field  data  in  tabular  form  for  easy  generation  of  3-D 
mesh  plots.  Zeros  are  returned  if  the  'elliptie'  routine  is  not  ealled.  The  format  is  as  follows: 


'earray.txf 

Rows  eorrespond  to  node  numbers  i,  while  eolumns  eorrespond  to  time  inerements. 
Column  1 :  x  eoordinate  of  eaeh  node  i,  at  t 

Column  2:  eleetrie  field  E.  [kV/em]  at  eaeh  node  i,  at  t  =  0 
Column  3:  eleetrie  field  E.  [kV/em]  at  eaeh  node  i,  at  t 

Column  n,+2\  eleetrie  field  E.  [kV/em]  at  eaeh  node  i,  at  t 


The  file  'earray.txf  eontains  the  vaeaney  eoneentration  data  in  tabular  form  for  easy  generation  of 
3-D  mesh  plots.  Zeros  are  returned  if  the  'parabolie'  routine  is  not  ealled.  The  format  is  listed  as 
follows: 
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’carray.txt’ 

Rows  correspond  to  node  numbers  i,  while  eolumns  eorrespond  to  time  inerements. 
Column  1 :  x  coordinate  of  eaeh  node  i,  at  t 

Column  2:  eoneentration  c-  [ppm]  at  eaeh  node  i,  at  t  =  0 
Column  3:  eoneentration  c.  [ppm]  at  eaeh  node  i,  at  t 

Column  n,+2\  eoneentration  [ppm]  at  eaeh  node  i,  at  t 


5.  Sample  Problem-BST  Film 


An  example  problem  is  diseussed  here  to  demonstrate  the  software  eapabilities  and  input  file 
format.  The  physieal  system  analyzed  here  is  a  thin  film  of  uniform  thiekness  7  =  100  nm 
(figure  1).  Isothermal  eonditions  are  assumed  with  6  =  298  K ,  a  temperature  at  whieh  undoped 
BST  remains  eubie  in  phase  for  molar  eoneentrations  of  Sr  greater  than  0.3  (Tinte  et  ah,  2004). 
Doping  with  Mg  further  lowers  the  Curie  point  (Cole  et  ah,  2003). 

Requisite  material  properties  are  listed  in  table  1 ,  deemed  representative  of  partieular 
eomposition  Bao.eSro^TiOs.  The  dieleetrie  eonstant  is  ehosen  as  representative  of  the  BST 
film  systems  of  interest  (Cole  et  ah,  2003),  while  thermodynamie  properties  c  and  fj  of  the  bulk 
substanee  are  used  (Todd  and  Lorenson,  1952).  The  vaeaney  diffusivity  d  is  eomputed  from 

D  =  JikgO  I  ez  =  dkj^O  /  ,  (65) 

where  D  is  a  thermally-aetivated  diffusivity  assoeiated  with  Fiek's  first  law,  ]u  is  the  drift 
mobility  of  ionie  eharges,  and  is  a  referenee  defeet  eoneentration.  We  assume  here,  from 

drift  mobility  data  on  100  nm-thiek  BST  films  (Zafar  et  ah,  1998),  a  eonstant  value  of 

=  2(10)  V“'m“'s“' ,  leading  to  the  value  of  d  reported  in  table  1.  Parameter  a  is  estimated 

from  eonsideration  of  the  erystal  strueture  and  lattiee  parameters,  while  the  value  of  z  indieates 
that  eaeh  oxygen  vaeaney  eontributes  a  eharge  of  magnitude  of  two  free  eleetrons. 


Table  1.  Properties  of  BST  film  at  298  K. 


Parameter 

Value 

Parameter 

Value 

500 

d 

6.24(10)'°  r' m-'  s"' 

C 

100  J  mol  '  K  ' 

a 

ia-29  3 

10  m 

fj 

115  J  mol  '  K^' 

z 

2.0 
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The  particular  initial-boundary  value  problem  includes  the  following  boundary  conditions, 
corresponding  to  shorted  or  ground  electrodes  at  the  film  boundaries,  which  themselves  are 
constrained  to  remain  impenetrable  to  vacancy  flow: 

^(0,t)  =  0,  ^(r,t)  =  0,  (66) 

e(0,t)  =  0,  Q{T,t)  =  Q.  (67) 

Initial  conditions  are  applied  as 

^(x,0)  =  10-^  (68) 

and 

p(x,0)  =  32043  C/m'.  (69) 


The  spatially-constant  value  of  j(x,0)  =  10^®  corresponds  to  an  initial  concentration  of 
Cq  =  1  ppm .  Additionally,  vacancy  conservation  constraint  (59)  is  activated  throughout  the 
analysis.  The  spatial  domain  is  decomposed  into  300  segments  of  equal  length  Ax  =  1/3  nm , 
with  a  total  of  301  nodes.  The  time  domain  is  chosen  as  0  <  t  <  1  //s .  Input  files  used  for  the 
computation  are  as  follows: 

'parameters. inp' 

mixed 

301 

l.Od-6 

l.Od-12 

100000 

l.Od-7 

298 

500 

0.0 

6.24dl0 

6.24dl0 

1.0d23 

l.Od-6 

2.0 

115.0 

100.0 

'electrostatic,  inp' 

potential 

0.0 

potential 

0.0 

initial  charge  distribution 
32043. 0~ 

32043.0 


32043.0 
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'diffusion.inp' 

flux 

0.0 

flux 

0.0 

yes 

1.0 

1.0 


1.0 

1.0 


The  solution  is  shown  in  figure  4,  with  concentration  and  electric  field  mesh  plots  generated 
from  output  files  'carray.txf  and  'earray.txt',  respectively.  Although  a  steady-state  concentration 
is  not  reached  in  the  short  duration  of  this  transient  simulation,  a  buildup  of  vacancies  at  the 
boundaries  x  =  0  and  x  =  T  =  100  nm  occurs,  as  the  charges  migrate  from  their  initially  uniform 
distribution  towards  the  boundaries  where  the  voltage  =  0  .  Vacancy  migration  has  a  negligible 

effect  on  the  local  electric  field  E,  however,  which  is  initially  linear  due  to  the  constant  initial 
charge  distribution  and  exhibits  an  average  value  of  zero. 


Figure  4.  Vacancy  concentration  and  electric  field:  Cq  -  1  Ppm  ,  V  -  0  . 

Vacancy  concentration  distributions  can  be  explained  upon  further  examination  of  chemical 
potential  //  of  equation  17.  The  term  pdy/jd^  establishes  an  energetic  penalty  for  large 
concentrations  and  forces  the  vacancies  to  diffuse  towards  a  uniform  state.  The  term 
accounts  for  electrostatic  interactions  due  to  charges  of  the  vacancies,  and  leads  to  their 
migration  towards  locations  where  the  potential  (j)  is  small.  The  final  distribution  of  vacancies  is 
thus  established  from  a  balance  of  the  contributions  of  these  two  terms  in  the  chemical  potential. 
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6.  Conclusions 


A  continuum  model  for  elastic  dielectric  semiconductors  with  mobile  charged  point  vaeancies 
has  been  presented.  The  complete  3-D  theoretical  framework  ineludes  the  physies  of  continuum 
elasticity,  electrostatics,  mass  diffusion,  and  charged  defeet  kinetics. 

A  restricted  version  of  the  theory  has  been  implemented  in  a  finite  differenee  code  allowing 
solution  of  the  elliptie  equations  of  electrostaties  of  dielectries  coupled  with  the  transient 
parabolie  equation  of  charged  diffusion.  The  numerieal  analysis  is  at  present  limited  to  one 
spatial  dimension. 

User  doeumentation  has  been  provided  for  the  computer  implementation.  Given  in  the  report  are 
deseriptions  of  the  source  code  strueture,  user  instructions,  and  representative  input  fdes  for 
analysis  of  BST  thin  films  containing  charged  oxygen  vaeancies.  The  source  eode  has  been 
ineluded  in  the  appendix. 
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Appendix.  Source  Code 
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Intentionally  left  blank. 


20 


program  semiconductor  ID 
c 

implicit  none 
c 

character ( len=50 )  problem  type 
c 

integer  i,j,k,l,m,n 

integer  nmax  !max  #  of  nodes  in  grid 

integer  nsteps  !max  *  of  timesteps 

integer  nout, kount, count  !#  steps  between  output  requests 

integer  ke,kc 
c 

real*8  t,tmax,dt 
real*8  pl,p2 
real*8  dl,d2 
real*8  z 
real *8  temp 
real*8  alpha 
real*8  beta 
real*8  L0,LC 
real*8  Ldot 
real*8  dx 
real*8  x 
real*8  pfree 
real*8  kb 
real*8  e 
real*8  s 
real*8  cp 
real*8  caveO 

real*  8  ,  pointer, dimension ( : ) 
real*  8 , pointer, dimension ( : ) 
real*  8 , pointer, dimension ( : ) 
real*  8 , pointer, dimension ( : ) 
real* 8 , pointer, dimension ( :  ) 
real* 8 , pointer, dimension (  :  , 
real* 8 , pointer, dimension (  :  , 
c 

500  format (12el4 . 5) 
c 

pfree  =  8 . 854187817d-12 
kb  =  1.3806505d-23 
e  =  1 . 60217653d-19 
c 

open (unit=l , f ile= ' parameters . inp ' , status= ' unknown ' ) 
read ( 1 , * ) problem_type 
read ( 1 , * ) nmax 
read ( 1 , * ) tmax 
read ( 1 , * ) dt 
read ( 1 , * ) nout 
read ( 1 , * ) LO 
read ( 1 , * ) temp 
read ( 1 , * ) pi 
read ( 1 , * ) p2 
read ( 1 , * ) dl 
read ( 1 , * ) d2 
read ( 1 , * ) alpha 
read ( 1 , * ) beta 


! output  flags 

! time, max  time,  delta-time 
! permittivity  constants 
!bulk  diffusivity  constants 
! valence  per  vacancy 
! temperature  (K) 

! convert  to  per  m^3  from  ppm 
! convert  to  vol.  fraction  from  ppm 
! initial  and  current  domain  length  (m) 

! velocity  of  domain  (m/s) 

Igrid  spacing  (m) 

! coordinate  (nm) 

! permittivity  of  free  space  (C/Vm) 

! Boltzmann's  constant  (J/K) 

! charge  of  an  electron  (C) 

! specific  entropy  (J/molK) 

! specific  heat  (J/molK) 

! initial  avg  concentration  (ppm) 

:dconc  ! cumulative  change  in  cone  (ppm) 
ideonst  ! relative  dielectric  constant 

lefield  lelectric  field  (V/m) 

ipotential  lelectric  potential  (V) 

:  charge  lelectric  charge  (C) 

I  : rearray  lelectric  field  output  matrix 

I  rrcarray  1  concentration  output  matrix 
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read ( 1 , * ) z 
read ( 1 , * ) s 
read ( 1 , * ) cp 
close  ( 1 ) 
c 

print*, 'Problem  problem  type 
print*, 'Domain  initial  length  ' , LO 
print*, 'Number  of  nodes  ' , nmax 
print*, 'Max  time  and  delta  t  ',tmax,dt 
print*, 'Bulk  diffusion  constants  ',dl,d2 
print*, 'Bulk  dielectric  constants  ',pl,p2 
print*, 'Valence  per  vacancy  ',z 
print* ,' Temperature  ' , temp 
print*, 'Entropy  ',s 
print* ,' Specif ic  heat ' , cp 

print*, 'Conversion  ppm  to  per  cubic  meter  ', alpha 
print*, 'Conversion  ppm  to  volume  fraction  ',beta 
c 

nsteps=  tmax/dt 
nsteps=  nsteps/nout+1 

print*, 'Number  of  output  steps  ',nsteps 
c 

allocate (dconc (nmax) ) 
allocate (dconst (nmax)  ) 
allocate (efield (nmax) ) 
allocate (potential (nmax) ) 
allocate (charge (nmax) ) 
allocate (earray (nmax, nsteps)  ) 
allocate (earray (nmax, nsteps) ) 
c 

t=0.0 

do  i=l,nmax 
dconc ( i ) =0 . 0 
dconst ( i ) =pl 
efield ( i ) =0 . 0 
potential ( i ) =0 . 0 
charge ( i ) =0 . 0 
end  do 
do  i=l,nmax 
do  j=l, nsteps 
earray ( i , j ) =0 . 0 
earray ( i , j ) =0 . 0 
end  do 
end  do 

Ldot=0 . 0 
LC=L0 

kount=nout 

count=nout 

ke=0 

kc=0 

c 

100  if (t . le . tmax)  then 
c 

dx=LC/ (nmax-1 . ) 
c 
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if (problem  type . eq . ' maxwell ' ) then 

call  elliptic (dconc, dx, temp, pi , p2 , z , e, charge, dt , pf ree, 

1  alpha, t, dconst, nmax, e fie Id, potential, 

2  nout, count, tmax, earray, nsteps, ke) 
end  if 

c 

if (problem  type . eq fick ') then 

call  parabolic (dconc, dconst, dl, d2, temp, nmax, dx, Ldot, s, cp, 

1  beta, t, dt, z, e, kb, pfree, alpha, efield, charge, potential, tmax, 

2  nout, kount, caveO, earray, nsteps, kc) 
end  if 

c 

if (problem_type . eq . ' mixed ' ) then 

call  elliptic (dconc, dx, temp, pi , p2 , z , e, charge, dt , pfree, 

1  alpha, t, dconst, nmax, efield, potential, 

2  nout, count, tmax, earray, nsteps, ke) 
c 

call  parabolic (dconc, dconst, dl, d2, temp, nmax, dx, Ldot, s, cp, 

1  beta, t, dt, z, e, kb, pfree, alpha, efield, charge, potential, tmax, 

2  nout, kount, caveO, earray, nsteps, kc) 
end  if 

c 

t  =  t  +  dt 
LC  =  LC+(Ldot*dt) 
go  to  100 
end  if 
c 

open (unit=5 , f ile= ' earray . txt ' , status= ' unknown ' ) 
open (unit=6, file= ' earray. txt ' , status= ' unknown ' ) 
x=0.0 

do  i=l,nmax 

write(5,500)x, (earray (i,j) ,j=l, nsteps) 
write(6,500)x, (earray (i,j) ,j=l, nsteps) 
x=x+dx* 1 . 0d9 
end  do 
close (5) 
close (6) 
c 

deallocate (dconc) 
deallocate (dconst) 
deallocate (efield) 
deallocate (potential) 
deallocate (charge) 
deallocate (earray) 
deallocate (earray) 
c 

end  program  semiconductor_lD 
c 

c################################################################### 

c 

subroutine  elliptic (dconc, dx, temp, pi , p2 , z , e, charge, dt , 

1  pfree, alpha, t, dconst, nmax, efield, potential, 

2  nout, count, tmax, earray, nsteps, ke) 
c 

implicit  none 
c 

character (len=50) bcleft, beright, dummychar 
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200 

c 


c 


c 


c 


c 


c 


c 


integer  i,j,k,l,m,n 

integer  nmax, adim, nout, count, nsteps,  ke 
integer, pointer, dimension ( : ) : : index 

real* 8  dx, temp, pl,p2, z,e,pfree, alpha, t , dt , tmax, tout 

real*8  valleft, valright 

real*8  dconc (nmax) 

real*8  dconst (nmax) 

real*8  chargeO (nmax) 

real*8  charge (nmax) 

real*8  potential (nmax) 

real*8  efield(nmax) 

real*8  x (nmax) 

real*8  earray (nmax, nsteps ) 

real*8  factor , energy, d 

real*  8 , pointer, dimension (:,:):: amatrix 
real* 8, pointer, dimension ( : ) : rbvector 
real*  8 , pointer, dimension ( : )  : : xvector 
format (16, 5el4 . 5) 

open (unit=2 , f ile= 'electrostatic. inp ' , status= ' unknown ' ) 

read (2 , * ) bcleft 

read (2 , * ) valleft 

read (2 , * ) be right 

read (2 , * ) valright 

read ( 2 , * ) dummychar 

do  i=l,nmax 

read (2, *) chargeO  (i) 

end  do 
close  (2 ) 

adim=nmax 

if (bcleft . eg . ' potential ' ) then 
adim=adim-l 

end  if 

if (beright . eq . ' potential ' ) then 
adim=adim-l 

end  if 

allocate (amatrix (adim, adim) ) 
allocate (bvector (adim) ) 
allocate (xvector (adim) ) 
allocate (index (adim) ) 

do  i=l,nmax 

charge (i) =charge0 (i) +e*z*dconc (i) *alpha 
X (i) =dfloat (i-1) *dx 
end  do 

do  i=l , adim 
do  j=l,adim 
amatrix ( i , j ) =0 . 0 
end  do 

amatrix ( i , i ) =-2 . 0 

k=i  +  l 

m=i-l 
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if ( i . It . adim) amatrix ( i , k) =1 . 0 
if ( i . gt . 1 ) amatrix ( i , m) =1 . 0 
end  do 
c 

f actor=-dx*dx/pf ree 
do  i=l , adim 

j=i 

if (be left. eq.  ' potential ' ) j  = j  +1 
bvector (i) =charge ( j ) *factor/dconst ( j ) 
end  do 
c 

if (be left. eq. ' potential ' )  bveetor ( 1 ) =bveetor ( 1 ) -vallef t 
e 

if (beright . eq . ' potential ' ) then 

bveetor (adim) =bveetor (adim) -valright 

end  if 
e 

if (belef t . eq . ' f lux ' )  then 
bveetor ( 1 ) =bveetor ( 1 ) -2 . *vallef t*dx 
amatrix ( 1 , 2 ) =2 . 0 
end  if 
e 

if (beright . eq .' flux ' )  then 

bveetor (adim) =bveetor (adim) +2 . *valright*dx 
i=adim-l 

amatrix (adim, i) =2 . 0 
end  if 
e 

eall  LU  Deeomp (adim, amatrix, index) 
eall  LU  BaekSub (adim, amatrix, index, bveetor) 
e 

do  i=l , adim 

j=i 

if (belef t . eq . ' potential ' ) then 
potential ( 1 ) =vallef t 
j=i  +  l 
end  if 

potential ( j ) =bveetor (i) 

if (beright . eq . ' potential ' )  potential (nmax) =val right 
end  do 
e 

k=nmax-l 
do  i=2 , k 
j=i-l 
m=i  +  l 

efield (i) = (potential ( j ) -potential (m) ) / (2 . *dx) 
end  do 

efield(l)= (potential ( 1 ) -potential (2 ) ) /dx 
efield (nmax) = (potential ( k) -potential (nmax) ) / dx 
if (belef t . eq . ' flux ') efield ( 1 ) =vallef t 
if (beright . eq . ' flux ' ) efield (nmax) =val right 
e 

energy=0 . 0 
do  i=l , k 

energy=energy+deonst (i) *efield (i) *efield (i) *dx 
end  do 
e 
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d=energy*dx*df loat (nmax-1) / (potential (1) -potential (nmax) ) **2 . 
c 

do  i=l,ninax 

dconst (i) =pl/ ( (l.+p2*efield(i) *efield(i) ) ** (0.3333333333)  ) 
end  do 
c 

tout=tmax-dt 

c 

if (count . eq . nout . or . t . gt . tout) then 

ke=ke+l 

count=0 

if (t . It . dt) then 

open (unit=3 , f ile= ' outelec . txt '  , 

1  status= ' unknown ' ) 

else 

open (unit=3 , f ile= ' outelec . txt ' , position= ' append '  , 

1  status= ' unknown ' ) 

end  if 
do  i=l,nmax 

write(3,200)i,t,x(i) , potential (i) ,efield(i) , charge (i) 
ear ray (i,ke)=efield(i) /l.OdS 
end  do 
close (3) 

print* Effective  dielectric  constant  ',d 
end  if 
c 

count=count+l 

c 

deallocate (amatrix) 
deallocate (bvector) 
deallocate (xvector) 
deallocate (index) 
c 

return 

end 

c 

c-################################################################### 

c 

subroutine  parabolic (dconc, dconst, dl, d2, temp, nmax, dx, Ldot, s, cp, 

1  beta, t, dt, z, e, kb, pfree, alpha, efield, charge, potential, tmax, 

2  nout, kount, caveO, carray, nsteps, kc) 
c 

implicit  none 
c 

character (len=50) bcleft, bcright, constraint 

integer  i,j,k,l,m,n 

integer  nmax, nout, kount, nsteps, kc 

real*8  dconc (nmax) 

real*8  dconst (nmax) 

real*8  efield (nmax) , charge (nmax) , potential (nmax) 

real* 8  dl , d2 , temp, dx, t , dt , z , e, kb, pfree, beta 

real*8  alpha, Ldot, s, cp, tmax, tout 

real*8  valleft, valright 

real*8  concO (nmax) 

real*8  cone (nmax) 

real*8  cdot (nmax) , grade (nmax) 

real*8  diff 1 (nmax) , diff2 (nmax) , q (nmax) , x (nmax) 
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real*8  denom, phi , dphi 

real*8  c, cdotmax, cdottemp, cave, caveO, LO, LC 

real*8  test, scale 

real*8  carray (nmax, nsteps ) 

300  format (i6, 5el4 . 5) 
c 

open (unit=4 , f ile= ' diffusion . inp ' , status= ' unknown ' ) 
read ( 4 , * ) bclef t 
read ( 4 , * ) vallef t 
read ( 4 , * ) bcright 
read ( 4 , * ) val right 
read ( 4 , * ) constraint 
do  i=l,nmax 
read ( 4 , * ) concO  ( i ) 
end  do 
close  (4) 
c 

if ( t . eg . 0 . 0 ) then 
cave0=0 . 0 
L0=0.0 
do  1=1, nmax 

cave0=cave0+dx*conc0  (i) 

L0=L0+dx 
end  do 

cave0=cave0/L0 
end  if 
c 

test  =  1 . /6 . *dx*dx/ (dl*kb*temp) 

c  if (dt . gt . test) print* ,' Warning :  dt  may  be  unstable !', test 

c 

do  1=1, nmax 
X (i) =dfloat (i-1) *dx 
cone (i) =conc0 (i) +dconc (i) 
c=conc  ( i ) 

call  chem  potential (c, temp, kb, alpha, s, cp, phi, dphi) 
diffl (i) =-dl*dphi 
diff2 (i) =d2*e*z*efield (i) 
end  do 
c 

k=nmax-l 
do  1=2 , k 
m=i  +  l 
n=i-l 

grade  (i)  =  (diffl  (m)  -diffl  (n)  )  /  (2  .  *dx) 
q(i)=  grade ( i ) +dif f 2 ( i ) 
end  do 
c 

if (bclef t . eq . ' potential ' ) then 

grade (1)  =  (-3. *diffl (1) +4 . *diffl (2) -diffl (3) ) /  (2 . *dx) 
q(l)=  grade ( 1 ) +dif f 2 ( 1 ) 
end  if 
c 

if (bcright . eq . ' potential ' ) then 
i=nmax-l 
j  =nmax-2 

grade (nmax) = (3 . *diff 1 (nmax) -4 . *diff 1 (i) +diff 1 ( j ) ) / (2 . *dx) 
q(nmax)=  grade (nmax) +diff2 (nmax) 
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end  if 


c 

if (bclef t . eq . ' flux ' ) q ( 1 ) =vallef t 
if (bcright . eq . ' flux ' ) q (nmax) =val right 
c 

if (bclef t . eq . ' special ' ) then 
c=conc  ( 1 ) 

call  chem  potential (c, temp, kb, alpha, s, cp, phi, dphi) 
q ( 1 ) =vallef t/ ( 1 . -cone ( 1 ) *beta) * 

1  ( phi +po tent ial ( 1 ) * charge (1) +0.5*efield(l) *efield(l) 

2  *dconst (1) *pfree+ (1 . /beta-conc (1) ) *dphi/beta+e*z*potential (1) * 

3  ( 1 . -cone (nmax) *beta) ) 
end  if 

c 

if (bcright . eq special ' )  then 
c=conc (nmax) 

call  chem  potential (c, temp, kb, alpha, s, cp, phi, dphi) 

q (nmax) =valright/ ( 1 . -cone (nmax) *beta) * (phi+potential (nmax) * 

1  charge (nmax) +0 . 5*efield (nmax) *efield (nmax) *dconst (nmax) *pfree 

2  + (1 . /beta-conc (nmax) ) *dphi/beta+e*z*potential (nmax) * 

3  ( 1 . -cone (nmax) *beta) ) 
end  if 

c 

k=nmax-l 
do  i=2 , k 
m=i  +  l 
n=i-l 

cdot (i)=- (q(m) -q(n) ) / (2.*dx) 
end  do 
c 

cdot (1) =- (-3. *q(l) +4 . *q(2) -q(3) ) / (2 . *dx) 
i=nmax-l 
j  =nmax-2 

cdot (nmax) =- (3 . *q (nmax) -4 . *q ( i ) +q ( j ) ) / (2 . *dx) 
c 

if (be left. eq. ' potential ' ) cdot ( 1 ) =0 . 0 
if (bcright . eq . ' potential ' ) cdot (nmax) =0.0 
c 

do  i=l,nmax 
cdottemp=cdot (i) 
cone (i) =conc (i) +cdot (i) *dt 
if  (conc(i) . It . 1 . d-15 ) cdot ( i ) =0 . 0 
cone (i) =conc (i) -cdottemp*dt 
end  do 
c 

if  (constraint . eq .' yes ') then 
cave=0 . 0 
LC=0.0 
do  i=l,nmax 

dconc (i) =dconc (i) +cdot (i) *dt 
cave=cave+ (concO (i) +dconc (i) ) *dx 
LC=LC+dx 
end  do 

cave=cave/LC 
scale=cave0/cave 
do  i=l,nmax 

dconc (i) =dconc (i) * s cal e- cone 0 (i) * (1. -scale) 
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end  do 
else 

do  i=l,ninax 

dconc (i) =dconc (i) +cdot (i) *dt 
end  do 
end  if 
c 

cdotmax  =  0.0 
cave  =  0.0 
n=0 

do  i=l,nmax 

if (abs (cdot ( i ) ) . gt . cdotmax) then 

cdotmax=abs (cdot (i) ) 

n=i 

end  if 

cave  =  cave  +  conc(i) 
end  do 

cave  =  cave/df loat (nmax) 
c 

if (bcright . eg . ' special ' . or . bclef t . eq . ' special ' ) then 
Ldot=q ( 1 ) / ( 1 . -cone ( 1 ) *beta) +q (nmax) / ( 1 . -cone (nmax) *beta) 
else 

Ldot=0 . 0 
endif 
c 

tout  =  tmax-dt 
c 

if ( kount . eq . nout . or . t . gt . tout) then 

kount=0 

kc=kc+l 

if (t . it . dt)  then 

open (unit=4 , f ile= ' outdif f . txt '  , 

1  status= ' unknown ' ) 

else 

open (unit=4 , f ile= ' outdif f . txt ' , position= ' append '  , 

1  status= ' unknown ' ) 

end  if 

print* Time t, ' C  average cave 
print*, 'Max  cone,  rate',  cdotmax,  n 

if (constraint . eq . ' yes ') print* , ' Constrained  rate : ' , test 
if (bcright . eq . ' special ' . or . bclef t . eq . ' special ' ) then 
print* , ' Ldot ' , Ldot 
end  if 
do  i=l,nmax 

write(4,300)i,t,x(i) , cone (i) , cdot (i) , q (i) 
car ray (i, kc) =conc  (i) 
end  do 
close  (4) 
end  if 
c 

kount=kount+l 

c 

return 

end 

c 

c-#################################################################### 

c 
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subroutine  chem  potential (c, temp, kb, alpha, s, cp, phi, dphi) 
c 

implicit  none 
c 

real*8  c, temp, kb, alpha, cp, s, phi, dphi 
real*8  gtotal 
real*8  xv 
real*8  avagadro,r 
real*8  nt 
c 

avagadro  =  6 . 022 14 1 9 9447d23 
r  =  8.314472 
c 

nt  =  alpha*1.0d6 
XV  =  c*alpha/nt 
c 

gtotal  =  temp* (cp-s+r*log (xv) ) 
c 

phi  =  nt/avagadro*gtotal 
c  dphi  =  kb*temp*c 

dphi  =  kb*temp*log (c) 
c 

return 

end 

c 

c-#################################################################### 

c 

subroutine  LU_Decomp (n, a, index) 
c 

implicit  double  precision  (a-h,o-z) 
dimension  a(n,n),  index (n),  v(n) 
c 

tiny  =  1 . Oe-2  0 
c 

do  i  =  l,n 
a  max  =  0.0 
do  j  =  l,n 

a_max  =  max (a_max, abs (a (i, j ) ) ) 
end  do  ! j 

V ( i )  =  1.0  /  a  max 
end  do  ! i 
c 

do  j  =  l,n 
c 

do  i  =  1 ,  j -1 
sum  =  a ( i , j ) 
do  k  =  1 , i-1 

sum  =  sum  -  a(i,k)  *  a(k,j) 
end  do 

a ( i , j )  =  sum 
end  do 
c 

a_max  =  0.0 
do  i  =  j , n 
sum  =  a ( i , j ) 
do  k  =  1 ,  j -1 

sum  =  sum  -  a(i,k)  *  a(k,j) 
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end  do 

a ( i , j )  =  sum 
dummy  =  v(i)  *  abs(sum) 
if  (  dummy  .gt.  a  max  )  then 
imax  =  i 
a  max  =  dummy 
end  if 
end  do 
c 

if  (  j  .ne.  imax  )  then 

do  k  =  1 , n 

dummy  =  a (imax, k) 

a ( imax, k)  =  a ( j , k) 

a ( j , k )  =  dummy 

end  do 

v(imax)  =  v(j) 
end  if 

index ( j )  =  imax 
c 

if  (  a(j,j)  .eq.  0.0  )  a(j,j)  =  tiny 
if  (  j  .ne.  n  )  then 
dummy  =  1.0  /  a ( j , j ) 
do  i  =  j  +1 , n 
a(i,j)  =  a(i,j)  *  dummy 
end  do 
end  if 
end  do  ! j 
c 

return 

end 

c 

c-################################################################## 

c 

subroutine  LU  BackSub (n, a, index, b) 
c 

implicit  double  precision  (a-h,o-z) 
dimension  a(n,n),  index (n),  b(n) 
c 

ii  =  0 
c 

do  i  =  l,n 

m  =  index ( i ) 

sum  =  b (m) 

b  (m)  =  b  ( i ) 

if  (  ii  .ne.  0  )  then 

do  j  =  ii,i-l 

sum  =  sum  -  a(i,j)  *  b(j) 

end  do 

else  if  (  sum  .ne.  0.0  )  then 
ii  =  i 
end  if 
b(i)  =  sum 
end  do 
c 

do  i  =  n, 1 , -1 
sum  =  b ( i ) 

if  (  i  .It.  n  )  then 
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do  j  =  i+l,n 

sum  =  sum  -  a(i,j)  *  b(j) 
end  do 
end  if 

b ( i )  =  sum  /  a ( i , i ) 
end  do 
c 

return 

end 

c 

c============================== 
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